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1. INTRODUCTION 


This report gives a brief summary of the research accomplished under this 
grant and a more detailed report on the activities during the last reporting 
period, from July, 1979 to January, 1980. The major objective of this re- 
search was to explore the application of singular perturbation theory for 
deriving algorithms suitable for on-board, real-time computation of optimal 
aircraft trajectory control. Minimum time Intercept In three dimensions was 
selected as a pilot problem formulation, and data for an early version F-4 
aircraft was used to represent aircraft aerodynamic and propulsion character- 
istics. The major products of this research are the derivation of nonlinear, 
near-optimal feedback control laws that are readily implementable in an air- 
craft flight computer, performance results for a variety of initial conditions 
(including optimal intercept at short ranges), and a point-mass three dimensional 
simulation suitable for representing fighter aircraft under closed-loop optimal 
control. A new grant has been initiated to utilize the results of this re- 
search first in a piloted simulation at NASA Langley and later in an actual 
flight test at the NASA Dryden research center. 

2. SUMMARY OF RESULTS FOR 1978 

The details of the research accomplished during the first 18 months are 
documented in References [1] - [3]. At this stage it was felt that a major 
milestone had been reached and Che results were published in [4]. Reference 
[1] presents the problem formulation and singular perturbation analvsis 
of the position, energy and heading dynamics. Altitude and flight path 
angle dynamics were not optimized. Reference [2] demonstrates the feasibility 
of on-board computation of the optimal control solution. During this perioa 


the optimal control solution was implemented in a feedback form as a control 
subroutine in a computer simulation of the F-4 aircraft. Since the altitude 
and flight path angle dynamics were not optimized, a simple proportional 
controller was used to calculate the required lift in the vertical plane 
based on the optimal altitude derived from optimization of the instantaneous 
energy, heading and relative position dynamics. Performance results were 
obtained for a variety of intercept conditions and an estimate was obtained 
of the CPU time and storage requirements for the control algorithm. In 
addition, the analysis was extended to the case of intercept at short ranges. 

This exhibits solutions that do not contain a constant altitude cruise leg. 

3. SUMMARY OF RESULTS FOR 1979 

The research this past year has been concentrated in four lin areas: 

1. Optimization of altitude and flight path angle dynamics 

2. Optimization for the short range case 

3. Implementation and evaluation of a proportional-plus-integral 
(P-I) controller for vertical lift calculation when altitude and 
flight path angle dynamics are not optimized. 

4. Off-line, open loop optimization using a conjugate gradient 
numerical optimization program. 

This section will detail the analysis for optimization of altitude and flight 
path angle dynamics and then conclude with example numerical results for the 
four areas cited above. 

3.1 Optimization of Altitude and Flight Path Angle Dynamics 

The analysis in this section assumes familiarity with the analysis in 
[!,■ 'Old [2]. The same notation is continued here. At this point the performance 
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index was generalized from minimum time co Include a term oenallzlne lar?e 
flight path angles 


- f 


1 + k sin Y dy 


( 1 ) 


and in addition a constraint was imposed on the maximum value for y 


1^1 1 Y_.„ < it/ 2 
max 


( 2 ) 


The above modifications to the problem formulation do not alter the analysis 
in [1] and [2]. 

Third Boundary Layer 

The Hamiltonian in the third (altitude) boundarv laver Is: 


■ W cos Y 


(3) 


H, - H, (h.E.S.Y) - X_ L KV/qsW + Xo L g/WV cos Y + ), V sin Y 
j X n °2 ^ ^ 


+ V, (Y^ -Y^) + V, cos^ Y - L‘} - 0 

1 max 2 max n 


(4) 


where is the first boundary layer Hamiltonian evaluated at present values 
of altitude (h) , energy (E), heading (0) and flight path angle 

Hj^ ■ (X^ cos 3 + X sin S) V cos Y “ X^ cos y^. 

^ 0^0 '0 


+ Y^ (T - D ) V/W + 1 + k sin‘- 

O 


(5) 


and 


D ■ q C_ + KW cos v/os 
0 3 D 

o 


( 6 ) 


The parameters X , X , X and X are known from previous outer and boundarv 

X V E, s« 

0 0 1 2 

layer solutions. The controls are y and the .norizontal component of lift 


I 
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So ve !iave 


JHj/Sv - 0 (7) 

5H./3L - 2X^ L KWosV: + r/TyT cos Y - 2v,L - 0 (fl) 

^ n K n H ■ J ^ 


Note Chat v >■ 0 off Che constraint bound and L can be determined as a 

n 

2 2 2 2 

function of y from (7). On the constraint bound L ■ L - wcos y and v 

n max 

can be determined as a function of y. In either case, L and v can be 


reearded as functions of y and other known oarameters. Hence, we need only 


combine (7) with (4) to compute the ootlmum flight path angle (y^) 
corresponding costate (X^ ). To do this we follow the procedure used In 


the first and second boundary/ layer and write 


Yj - ar? (mjx ) . ,1^ (h^ - h} 

where Is Che second boundary layer Hamiltonian 
- H - X \} KV/osW + X- L g/W cos y 


(^) 


( 10 ) 


Tne maximization in (9) Is performed over the range 0 < y < v , where we 

“ nax 

have used the fact that sign {X. } ■ - sign (y) ■ sign {h - h«} and h. is 

h 2 1 

the optimum altitude from the second boundair/ laver. Since It is possible 
for to change sign it was decided to instead use the form 


- arg {min [H^/sin y]} • sign {h 2 - h} 


Having determined y. and L , >, can be computed using (4) 
3 P3 h3 

V * -[H, - L‘ KV/osW + X, L g/W cos v ]/’’ sin v 
I13 1 '^3 "3 ^ 


( 11 ) 


( 12 ) 
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Numerical results for as a function of the altitude error for a 
zero and non-zero heading errors are shown in Fig. 1. The effect of 
the weighting parameter (k) in the performance index (1) is also shovm. 


Fourth Boundary Laver 

For solutions off the lift constraint bound we write the fourth boundary 
laver Hamiltonian in the form 


H 


, - H.(h,E,3,Y) - K(25L,W cos y + 6LF + 2(6L L 

4 3 c.. 1 X n r 


+ *L")V/asW, + ,\. £L g/U’V cos y + X 5L g/T-T - 0 (13) 

n ■ D <> n y 1 

where 5L, and £L are perturbations in the controls from the third boundarv 
In 

laver solutions. 


- '.•/ cos y 


oL ■ L - L 
n n n. 


(14) 

(15) 


Using the conditions H, ■ 0, 3H,/35L, ■ 0 and 3H,/3£L ■ 0 we obtain the 

4 4 1 4 n 

following solutions : 


5L 


n, 


gns/Z'^ KV“ cos V - L 


-:l - r-H qsU/X., ]'■' • sign (y, - v) 

1 4 


(16) 

(IT) 


The final lift and bank angle equations are 
L 


\ 

(V cos ■' £L, ,)" + L*” ]■ 

1-* n, 

4 


■j - can ^(L /L , ) 
n, 14 


(13) 


(19) 


If L in (13) exceeds L Chen we set L ■ L and refomulate the necessarv 

max nax 

conditions to find Che optimum bank ancle. Lettlnc 


L, ■ L cos u 
1 max 


L > L sin u 
n max 


Then from 3H,/3u ■ 0 we obtain 

4 


( 20 ) 

( 21 ) 


can u ■ (1- /X cos y) • sign (* - ♦) 

02 o 


Using (22) in K, - 0 we can evaluate X as one of the roots of 
4 Y 


AX" + BX + C - 0 
V Y 


where 


o 2 2 2 

A ■ g""[L /!v - cos" y]/V 

max 


B ■ 2® cos yR/V 


C - (X L g/T.-V cos. Y)^ - 
^2 >nax 


( 22 ) 


(23) 


(24) 

(25) 


(26) 


and 

5 - H - X_ (L^ - W^cos\) KV/osM + X. ” sin •' (27) 

X ti^ riAX ii ^ 

In Appendix A it is shown that (23) will always have real roots with opposite sign. 
So the root of (23) is chosen such that 


sipn {X_^} - -sign - y} 


(2S) 


'.’ith this value of the ootimuo bank ancle is comnuted using (22) 
care to place u in the proper quadrant based ^.so on Che sign of (y 


3 


takina 

- y) 


3.2 Numerical Results 


This section documents the results of numerical testa made to verify the 
analytical results In this report. The aircraft simulation model used in this 
study Is the same as In [2]. The optimal control computation requirements. 
Including the third and fourth boundary layer optimization, were 0.07 seconds 
for update and 50198 bytes of core space versus 0.03 seconds per update and 
43724 bytes for the suboptlmal control described In [2]. The Increase In com- 
putation time l.« primarily due to the third boundary layer computation, which 
requires a numerical search over (11). 

Two test cases were Investigated In the development of the fully optimized 
control. The first test :ase In Case I from [2] starting at t ■ 87 seconds, 
at which point heading error Is near zero, with the aircraft moved to an altitude 
of 40,000 leer while holding total Initial energy constant, thus Introducing 
a large Initial altitude error. This case Is hereafter referred to as Case I-A. 
The second test case examined Is Case III In [2]. Case III was chosen for the 
hard Initial turn required of the aircraft for Intercept, which resulted In a 
commanded altitude above the climb path. The aircraft In this case Is Initially 
close to the zero-heading-error optimum climb path. The Initial geometries In 
the horizontal plane are displayed In Fig. 2. The trajectories were stopped at 
200 seconds for Case I-A and at 60 seconds for Case III, displaying Initial trans- 
ient phases only. This Is sufficient to show the effect < the third and fourth 
boundary layer optimization since the control solution converges to the second 
boundary layer solution as the flight path angle ind altitude errors approach 
zero. Figures 3 through 15 display cctunl and commanded flight path angle and 
altitude histories and lift and bank angle histories for Cases I-A and III and 
altitude versus aach number plots for Case III, both with fully optimized control, 
and with proportional vertical lift control as in [2]. 


n' 


1 

I In Case I-A the control failed twice due to numerical difficulties encoun- 

tered when the altitude error became small. This can be seen in Fig. 7 where 
I the commanded flight path angle (GMD) discontinuously jumps to its constrained 

I maximum at approximitely 90 and 160 seconds. Note the altitude errors in Fig. 4 

at these times. The computer implementation of the control algorithm called 
j for this when the first bracketed term of (11) became negative. The numerical 

problem leading to the excursions can be traced to the fact that both the num- 
erator and denominator in (11) asymptotically approach aero as the aircraft al- 
titude approaches the second boundary layer commanded altitude (the optimal 
flight path angle also approaches zero under these circumstances) . It should 

I be noted that the calculation of in [2] takes the form: 

« 2 

X. - (X_ V-KL^^/qsW - mV./L (29) 

I fc n X ^ ti 

^ and is prone to considerable inaccuracy when the heading error is very small, 

] since approaches 0. Doubtless, the carryover of second boundary layer in- 

accuracies into the thirc. aggravated the numerical difficulties. It can be 
^ shown, however, employing I'Hospital's rule, that (11) does in fact converge 

to zero. W- this encouragement, modifications are currently being made to 
■" the algorithm to circumvent these problems. 


No discontinuities are seen in the Case III trajectory (Figs. 9, 11, 13) 
since the third boundary layer control solution had not converged to that of 
the second at any point. It is of interest in this case to note the appreciable 
difference between the altitude versus mach plots for proportional and fully 
optimized control, seen in Figs. 14 and 15. 

Although no test trajectories were run fully to intercept, it is estimated, 
based on use of the tabular first boundary layer data in [2] over straight-line 
approximations for the intercept paths, that the time savings accrued in both 


- 9 - 


TGURE 2. 


TARGET 



AIRCRAFT 


CASE I-A 


TARGET 

•4 


1 AIRCRAFT 

4 ► 


CASE III 


HORIZONTAL PLANE PROJECTIONS OF INITIAL GEOMETRIES FOR CASES 
TREATED IN STUDY 


- 10 - 






13.00 ^ 0.00 80.00 120.00 160.00 200.00 240.00 

TME (SEG 

FIGURE 4. ACTUAL AND COMMANDED ALTITUDE FOR CASE I-A WITH FULLY OPTIMIZED 


LIFT AND BANK ANGLE 







^.00 40.00 80.00 120.00 160.00 200.00 240.00 

TME (SEC] 

FIGURE 5. LIFT AND BANK ANGLE FOR CASE I-AWITH FULLY OPTIMIZED LIFT AND 


BANK ANGLE 



0.00 


> 40.00 


160.00 


200.00 


2 ^ 0.00 


FIGURE 6. 


80.00 120.00 

TME (SEC) 

ACTUAL AND COMMANDED FLIGHT PATH ANGLE FOR CASE I-A WITH 


SUBOPTIMAL PROPORTIONAL CONTROL 






GMD (DEG) 












00*08 













FIGURE 14. ALTITUDE VS MACH ^^L’MBER FOR CASE III WITH SUBOPTLMAL PROPORTIONAL CONTROL 



P 




Cases for full control optimization Is on the order of 2 seconds. An exact 
comparison will be made after modifications have been completed to avoid num- 
erical difficulties, as cited above. 

The development In [2] for range-matching to achieve optimal short range 
Intercept was tested with initial geometries similar to Cases I and III from 
[2], placing the target at several values of Initial range. The target altitude 
was set at 20,000 feet for these tests. Horizontal plane plots of the resulting 
trajectories are displayed In Figs. 16 to 21. The altitude histories and alti- 
tude versus mach plots for these tests are exhibited In Figs. 22 to 26. 

An attempt was made to Improve the aircraft's suboptlmal tracking of optimal 
altitude commands through the Implementation of a proportlonal-plus-lntegral 
scheme for vertical lift. This was motivated by the Inability of the suboptlmal 
proportional controller to accurately follow the first boundary layer descent 
path. For this aircraft, the optimal descent path lies along the maximum mach 
number boundary illustrated by the dashed line In Fig. 24. For the suboptlmal 
control, this corresponds to a ramp In the commanded altitude requiring a type 
one controller for nulling the steady-state error. It Is obvious from this 
figure that a proportional controller results in a steady-state altitude error 
during descent. 

Initially, a simple addition of an integrator on altitude and error was 
added to the proportional scheme already In use In the manner described In 
Appendix B. This controller was tested on Case I-A. Profiles of actual and 
commanded altitude as functions of time are shown In Fig. 27 for the original 
proportional controller, while Fig. 23 reflects the inclusion of the integrator. 

The considerable instability exhibited in Fig. 28 was due to coupling between 
flight path angle (y) and energy rate (E) dynamics described more fully in .Appendix 
B. On the basis of this problem, the controller was reformulated with time 


constants compensated as functions of energy in the manner described In Appendix 
B. This formulation was tested In the same manner as the first, resulting In 
the altitude profile displayed In Fig. 29. With the encouragement of enhanced 
stability and greatly reduced errors, the controller was tested on Case I-A to 
Intercept. Unfortunately, the controller became completely unstable at cruise 
due to large coupling between energy rate dynamics and lift In that region, 
this being due to an extremely low energy rate. Weighing the difficulties antic- 
ipated In solving this problem against the expected benefits, no further work 


will be done In this area. 
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APPENDIX A 


This aroendix oroves that the roots of the Quadratic exnression for X 

... y 

obtained in (23) of section 3.1 are alwavs real and of onnosite sien. It 
will also be shown that the fourth boundary layer controls ara asymptotic 
to those of the third. 

From section 3.1 we have that, on the lift constraint 


AX^ + BX^ + C - 0 (A-1) 

Y Y 

where 

A = 3“[L ~nr - cos^ y]/v" (A-2) 

max 

3 = 2$g cos Y/V (A-3) 

C = [X L g/n: cos y]^ - (A-4) 

a., max 

? = H, (h.E.S.Y.L =■ L ) + X. V sin V (A-5) 

1 max h^ 

where has been defined in (5). 

The roots of (A-1) are as follows: 

X = [-B + /b" - 4 AC ]/2A (A-6) 

^ 1,2 

where the term under the radical is 


B~ - 4aC 


4[('JL - 

max 


(-5L s/n:) 

max' 


+ (X„ (L cos y)^] 

pj ^ TTlaX 


O 70 

(X, L )“ 

c, max' 


CA-7) 


For the value of this expression to remain non-neeative it is necessarv 



I ^ 


* 


r + [1 - (L /W cos y)“] [\„ <i/v]“ > n 
nax ^ a ‘ — 


(A-P) 


or 


J" > [L “ - T" cos" y] [X- cos y]‘ 

— tWQV R “ 


max 


(A-Q) 


In Che Chird boundarv laver it can be seen that 


2 2 2 2 
L - L - W cos Y 

max 


(A-10) 


Also, the Chird boundary layer Hamiltonian in (4), evaluated at the con- 
ditions under consideration is 


K ■ + Xg L cos V ^ 0 

^ “2 max 


(A-11) 


The product X. . L is always negative, so (A-10) can be rewritten 
32 ^ 


<f > I X„ L gA’T-T cos y| 

— ' 6~ n ‘ ' 

^ max 


(A-12) 


Combining (A-10) and (A-12) gives us (A-9) assuring us that the roots of (A-1) 
are real. 

Adressing the signs of the roots (A-6) , it is obvious that, in the 
quadratic formula, the term under the radical must have a dominating magnitude 
for two different signs to occur. This is tantamount to requiring Chat 
AC < 0 (A-13) 


Since it can be seen from (A-2) that A > 0 where L > I’, then we must have that 

— max — 


C <_ 0 (A-14) 

From (A-4) and (A-12) it can be seen that this is indeed the case. 

Finallv, to prove the as''rmntotic behavior of the fourth boundarv laver 
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controls, we will show that oLi and 5L both tend to zero as the current fllcht 

n 

path angle approaches its third boundair/ layer optimal value fTori the 

fourth boundary laver analysis we have that, off the lift constraint 




[X qsg/2X KV“ cos Y - L_ 1 

El n^ 


[-H,qsW/X„ KV - 6~L ]'" 

3 \ ^4 


(A-15) 


(A-16) 


We note from (8) that values of L satisfying 3H/3L ■ 0 off the lift 

n^ n 

constraint are 

L -X csg/2X KV“ cos y (A-17) 

n. 3, E 3 


It is evident from the above that 6L in (A-15), being continuous, will 
asNTnpotically approach zero as the current flight path angle (y) v^. 

Turning to in (A-16) it seen that 5L^ and will both approach 


zero y 


APPENDIX B 


This appendix documents the derivation of a proportional -nlus-lntegral 

scheme for sub-optimally computing vertical lift in the simulation. Figure B-1 

illustrates a straightforward addition of a term integrating (h - h) to the 

c 

proportional controller shown in Figure 6 of [2]. This approach proved un- 
satisfactory due to the interaction of vertical lift with aircraft energy 
rate, which was not accounted for, and which resulted in unstable behavior. 

The dashed line in Figure 1 models the effect of vertical lift 

energy' rate (E) , and is seen to act as an outer loop around the original 

model for the third order P-I controller. 


On 

the basis of this the control gains K^., and K^, 

defined by"*" 

K, 

” 1/ T_ 


I 

I 

(B-1) 

K 

- 1/t 

(B-2) 

Y 

Y 




(B-3) 


were redisigned to account for E dvnamics. Tlie term h^(E) represents the 
optimum climb (descent) altitude schedule as a function of energy level. 

A linearized perturbation block diagram is shown in Figure B-2. The 
linearization is about L » V., where 


3h /3E 

(B-4) 

c 

-3E/3E - 2KV/qs 

(B-5) 


The characteristic equation for this system is given by 


and here correspond to and in Figure 6 of [2]. 
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(B-6) 


3 -> 

s + [K^ + aK^K^ls + ''’ ° 


where 


-n^K| 


cos Y 


(B-7) 


If we arrange Che poles of the system according Co 


(s + K^) (s + KjS + K 3 ) - 0 


(B-8) 


where 


K, > 10 u 

1 — n 


(B-9) 


Chen we can compute Che gains equating Che coefficients of like powers 
of s in (B-6) and (B-8) : 


K, - Ki + K2 - a(K^K2 + K3) + 


\ “ ^"iV "3 - ""l"3>^ 


= K K /K 
I 1 3 Y 


Note that "a" from (B-7) is undated continuously durini? a trajectorv and 
renresents a compensation tern to account for E d'mamics. Setting a ■ 0 


reduces the solution to the form illustrated in Figure B-1. 
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